﻿using GeometryCL;
using GeometryCL.GeoAPI.Geometries;
using NLog;
using NTSDecimalCL;
using NTSDecimalCL.NetTopologySuite.Algorithm;
using NTSDecimalCL.NetTopologySuite.Triangulate.QuadEdge;
using ProtoBuf;
using System;
using System.Collections.Generic;
using System.IO;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace NTSDecimalCL.DivideAndConquerALG
{
    class DivideAndConquerAlg
    {
        private QuadEdgeSubdivision _sv;
        private ICollection<Vertex> S;

        private static NLog.Logger logger = LogManager.GetCurrentClassLogger();

        public DivideAndConquerAlg(QuadEdgeSubdivision sv, ICollection<Vertex> S)
        {
            _sv = sv;
            this.S = S;
        }

        public void Execute()
        {
            var frame = _sv.GetFrameVertex();
            foreach (var Vertex in frame)
            {
                S.Add(Vertex);
                //logger.Info("frame Vertex:" + Vertex.ToString());
            }
            _sv.DeleteAllEdge();
            S = S.OrderBy(vertex => vertex.X).ThenBy(vertex => vertex.Y).ToList();
            S = RemoveTheSame(S);
            DelaunayAlgReturn r = Delaunay(S, 1);
            _sv.SetStartingEdge(r.le);
            //TimeComputeCL.TimeCompute.Log();
        }

        private ICollection<Vertex> RemoveTheSame(ICollection<Vertex> S)
        {
            ICollection<Vertex> r = new List<Vertex>();
            Vertex last = null;
            foreach(var vertex in S)
            {
                if (last == null || !vertex.Equals(last))
                {
                    r.Add(vertex);
                }
                last = vertex;
            }
            return r;
        }

        private DelaunayAlgReturn Delaunay(ICollection<Vertex> S, int h)
        {
            //if (S.Count == 11)
            //    logger.Info("debug");            
            //logger.Info("S.Count:" + S.Count + ";h:" + h);
            try
            {                
                if (S.Count == 2)
                {
                    QuadEdge a = _sv.MakeEdge(S.ElementAt(0), S.ElementAt(1));
                    //if (a.Orig.X == 7560 && a.Orig.Y == -7560 || a.Dest.X == 7560 && a.Dest.Y == -7560)
                    //    logger.Info("S.Count == 2:a:" + a);
                    return new DelaunayAlgReturn(a, a.Sym);
                }
                else if (S.Count == 3)
                {
                    // Let s1, s2, s3 be the three sites, in sorted order
                    Vertex s1 = S.ElementAt(0);
                    Vertex s2 = S.ElementAt(1);
                    Vertex s3 = S.ElementAt(2);
                    //bool islog = false;
                    /*if (s1.X == 7560 && s1.Y == -7560 || s2.X == 7560 && s2.Y == -7560
                        || s3.X == 7560 && s3.Y == -7560)
                    {
                        logger.Info("S.Count == 3:s1:" + s1 + ";s2:" + s2 + ";s3:" + s3);
                        islog = true;
                    }*/
                    // Create edges a connecting s1 to s2 and b connecting s2 to s3
                    QuadEdge a = _sv.MakeEdge(s1, s2);
                    QuadEdge b = _sv.MakeEdge(s2, s3);
                    QuadEdge.Splice(a.Sym, b);
                    // Now close the triangle
                    if (IsCCW(s1, s2, s3))
                    {
                        QuadEdge c = _sv.Connect(b, a);
                        /*if (QuadEdge.isget)
                        {
                            logger.Info("get1");
                            QuadEdge.isget = false;
                        }*/
                        //if (islog)
                        //    logger.Info("S.Count == 3:out1");
                        return new DelaunayAlgReturn(a, b.Sym);
                    }
                    else if (IsCCW(s1, s3, s2))
                    {
                        QuadEdge c = _sv.Connect(b, a);
                        /*if (QuadEdge.isget)
                        {
                            logger.Info("get2");
                            QuadEdge.isget = false;
                        }*/
                        //if (islog)
                        //    logger.Info("S.Count == 3:out2");
                        return new DelaunayAlgReturn(c.Sym, c);
                    }
                    else
                    {
                        // The three points are collinear
                        //if (islog)
                        //    logger.Info("S.Count == 3:out3");
                        return new DelaunayAlgReturn(a, b.Sym);
                    }
                }
                else
                {
                    // |S| >= 4. Let L and R be the left and right halves of S
                    ICollection<Vertex> L, R;
                    SplitList(S, out L, out R);
                    DelaunayAlgReturn ld = Delaunay(L, h + 1);
                    DelaunayAlgReturn rd = Delaunay(R, h + 2);                                        
                    QuadEdge ldo = ld.le;
                    QuadEdge ldi = ld.re;
                    QuadEdge rdi = rd.le;
                    QuadEdge rdo = rd.re;
                    // Compute the lower common tangent of L and R
                    do
                    {
                        if (rdi.Orig.LeftOf(ldi))
                            ldi = ldi.LNext;
                        else if (ldi.Orig.RightOf(rdi))
                            rdi = rdi.RPrev;
                        else
                            break;                        
                    }
                    while (true);

                    // Create a first cross edge basel from rdi.Org to ldi.Org
                    QuadEdge basel = _sv.Connect(rdi.Sym, ldi);
                    /*if (QuadEdge.isget)
                    {
                        logger.Info("get3");
                        QuadEdge.isget = false;
                    }*/
                    if (ldi.Orig == ldo.Orig)
                        ldo = basel.Sym;
                    if (rdi.Orig == rdo.Orig)
                        rdo = basel;
                    do // This is the merge loop
                    {
                        //bool islog = false;
                        //{LINESTRING(112.15865 21.81836, 112.15814 40.23361)}
                        //{LINESTRING(112.15721 32.06601, 112.1572 32.21389)}
                        //if (basel.id == 713)
                        //    islog = true;
                        // Locate the first L point(lcand.Dest) to be encountered by the rising bubble
                        // and delete L edges out of basel.Dest that fail the circle test
                        QuadEdge lcand = basel.Sym.ONext;
                        //if (islog)
                        //{
                        //    basel.Sym.LogONextEdge();                            
                        //}
                        if (Valid(lcand, basel))
                        {
                            //while (lcand.ONext.Dest.IsInCircle(basel.Dest, basel.Orig, lcand.Dest))
                            while (IsInCircumcircle(
                                basel.Dest.Coordinate, basel.Orig.Coordinate, lcand.Dest.Coordinate, 
                                lcand.ONext.Dest.Coordinate))
                            {
                                /*Vertex A = basel.Dest;
                                Vertex B = basel.Orig;
                                Vertex C = lcand.Dest;
                                Vertex D = lcand.ONext.Dest;*/
                                QuadEdge t = lcand.ONext;
                                //if (lcand.Orig.X == 0 && lcand.Orig.Y == 7560 && lcand.Dest.X == 7560 && lcand.Dest.Y == -7560
                                //|| lcand.Dest.X == 0 && lcand.Dest.Y == 7560 && lcand.Orig.X == 7560 && lcand.Orig.Y == -7560)
                                //logger.Info("DeleteEdge1:" + lcand);
                                //string strlcand = lcand.ToString();
                                _sv.Delete(lcand);
                                /*if (QuadEdge.isget)
                                {
                                    logger.Info("get4:lcand:" + strlcand + ";lcand.ONext:" + t);
                                    QuadEdge.isget = false;
                                }*/
                                lcand = t;
                            //    if (islog)
                            //        logger.Info("lcand2:" + lcand + ";A:" + A + ";B:" + B
                            //            + ";C:" + C + ";D:" + D);
                                if (!ldo.IsLive || !rdo.IsLive)
                                {
                                    logger.Error("Delaunay IsInCircle1 error");// :A:" + A + ";B:" + B + ";C:" + C
                                        //+ ";D:" + D);
                                    throw new Exception("IsInCircle1");
                                }                                
                            }                            
                        }

                        // Symmetrically, locate the first R point to be hit, and delete R edge
                        QuadEdge rcand = basel.OPrev;
                        //if (islog)
                        //    basel.LogOPrevEdge();
                        if (Valid(rcand, basel))
                        {
                            //while (rcand.OPrev.Dest.IsInCircle(basel.Dest, basel.Orig, rcand.Dest))
                            while (IsInCircumcircle(
                                basel.Dest.Coordinate, basel.Orig.Coordinate, rcand.Dest.Coordinate,
                                rcand.OPrev.Dest.Coordinate))
                            {
                                /*Vertex A = basel.Dest;
                                Vertex B = basel.Orig;
                                Vertex C = rcand.Dest;
                                Vertex D = rcand.OPrev.Dest;*/
                                QuadEdge t = rcand.OPrev;
                                //if (lcand.Orig.X == 0 && lcand.Orig.Y == 7560 && lcand.Dest.X == 7560 && lcand.Dest.Y == -7560
                                    //|| lcand.Dest.X == 0 && lcand.Dest.Y == 7560 && lcand.Orig.X == 7560 && lcand.Orig.Y == -7560)
                                    //logger.Info("DeleteEdge2:" + lcand);
                                _sv.Delete(rcand);
                                /*if (QuadEdge.isget)
                                {
                                    logger.Info("get5");
                                    QuadEdge.isget = false;
                                }*/
                                rcand = t;
                            //    if (islog)
                            //        logger.Info("rcand2:" + rcand + ";A:" + A + ";B:" + B 
                            //            + ";C:" + C + ";D:" + D);
                                if (!ldo.IsLive || !rdo.IsLive)
                                {
                                    logger.Error("Delaunay IsInCircle2 error");//:A:" + A + ";B:" + B + ";C:" + C
                                        //+ ";D:" + D);
                                    throw new Exception("IsInCircle2");
                                }                                
                            }                            
                        }
                        /*if (h == 2)
                            logger.Info("lcand:" + lcand);
                        if (h == 2)
                            logger.Info("rcand:" + rcand);*/

                        // If both lcand and rcand are invalid, the basel is the upper common tangent
                        if (!Valid(lcand, basel) && !Valid(rcand, basel))
                        {                            
                            //CheckIfError(S, basel);
                            break;
                        }

                        // The next cross edge is to be connected to either lcand.Dest or rcand.Dest
                        // If both are valid, then choose the appropriate one using the InCircle test
                        if (!Valid(lcand, basel) || (Valid(rcand, basel) &&
                            //rcand.Dest.IsInCircle(lcand.Dest, lcand.Orig, rcand.Orig)
                            //NetTopologySuite.Triangulate.QuadEdge.TrianglePredicate.IsInCircleNonRobust(
                            IsInCircumcircle(
                                lcand.Dest.Coordinate, lcand.Orig.Coordinate, rcand.Orig.Coordinate,
                                rcand.Dest.Coordinate)
                            ))
                        {
                            // Add cross edge basel from rcand.Dest to basel.Dest
                            //QuadEdge baselSym = basel.Sym;
                            basel = _sv.Connect(rcand, basel.Sym);                            
                            /*if (QuadEdge.isget)
                            {
                                logger.Info("get6");
                                QuadEdge.isget = false;
                            }*/
                        }
                        else
                        {
                            // Add cross edge basel from basel.Org to lcand.Dest
                            //bool bget = QuadEdge.isget;
                            //QuadEdge baselSym = basel.Sym; 
                            basel = _sv.Connect(basel.Sym, lcand.Sym);                            
                            /*if (!bget && QuadEdge.isget)
                            {
                                logger.Info("get7:baselSym:" + baselSym + ";lcand.Sym:" + lcand.Sym
                                    + ";basel:" + basel);
                                QuadEdge.isget = false;
                            }*/
                        }
                        //if (h == 2)
                        //    logger.Info("basel:" + basel);
                    }
                    while (true);                    
                    return new DelaunayAlgReturn(ldo, rdo);
                }
            }
            catch(Exception e)
            {
                logger.Error("Delaunay error:" + e.Message + ";S.Count:" + S.Count);
                /*foreach(var vertex in S)
                {
                    logger.Error("X:" + vertex.X.ToString("G17") + ";Y:" + vertex.Y.ToString("G17"));
                }*/
                throw new Exception(e.Message);
            }
        }

        private bool IsCCW(Vertex p1, Vertex p2, Vertex p3)
        {
            return (p2.X - p1.X) * (p3.Y - p1.Y)
                 - (p2.Y - p1.Y) * (p3.X - p1.X) > 0;
        }

        public static void SplitList(ICollection<Vertex> S, out ICollection<Vertex> L, out ICollection<Vertex> R)
        {
            L = new List<Vertex>();
            R = new List<Vertex>();
            int mid = S.Count / 2;
            for(int i = 0; i < S.Count; i++)
            {
                if (i < mid)
                {
                    //if (S.ElementAt(i).X == 0 && S.ElementAt(i).Y == 7560)
                    //    logger.Info("SplitList left i:" + i + ";S.Count:" + S.Count);
                    L.Add(S.ElementAt(i));
                }
                else
                {
                    //if (S.ElementAt(i).X == 0 && S.ElementAt(i).Y == 7560)
                    //    logger.Info("SplitList right i:" + i + ";S.Count:" + S.Count);
                    R.Add(S.ElementAt(i));
                }
            }
        }

        private bool Valid(QuadEdge e, QuadEdge basel)
        {
            return e.Dest.RightOf(basel);            
            /*double ccw = (basel.Dest.X - e.Dest.X) * (basel.Orig.Y - e.Dest.Y)
                 - (basel.Dest.Y - e.Dest.Y) * (basel.Orig.X - e.Dest.X);
            if (ccw > 0)
                return true;
            else if (ccw < 0)
                return false;
            else
            {
                if (e.Orig.Y + e.Dest.Y > basel.Orig.Y + basel.Dest.Y)
                    return true;
                else
                    return false;
            }*/
        }

        private void CheckIfError(ICollection<Vertex> S, QuadEdge basel)
        {
            //TimeComputeCL.TimeCompute.Start("clist.Add");
            List<Coordinate> clist = new List<Coordinate>();
            foreach (var vertex in S)
            {
                clist.Add(vertex.Coordinate);
            }
            //TimeComputeCL.TimeCompute.Pause("clist.Add");
            //TimeComputeCL.TimeCompute.Start("clist.ToArray");
            NetTopologySuite.Geometries.GeometryFactory gf = new NetTopologySuite.Geometries.GeometryFactory();
            var arr = clist.ToArray();
            //TimeComputeCL.TimeCompute.Pause("clist.ToArray");
            //TimeComputeCL.TimeCompute.Start("new.ConvexHull");
            ConvexHull ch = new ConvexHull(arr, gf);
            //TimeComputeCL.TimeCompute.Pause("new.ConvexHull");
            //TimeComputeCL.TimeCompute.Start("get.ConvexHull");
            IGeometry g = ch.GetConvexHull();
            //TimeComputeCL.TimeCompute.Pause("get.ConvexHull");
            //TimeComputeCL.TimeCompute.Start("other");
            try
            {
                Coordinate[] pointList = g.Coordinates;
                if (pointList.Length > 1)
                {
                    bool isget = false;
                    for (int i = 1; i < pointList.Length; i++)
                    {
                        if (pointList[i - 1].Equals2D(basel.Orig.Coordinate) &&
                            pointList[i].Equals2D(basel.Dest.Coordinate) ||
                            pointList[i - 1].Equals2D(basel.Dest.Coordinate) &&
                            pointList[i].Equals2D(basel.Orig.Coordinate) ||
                            IsOnLine(basel.Orig.Coordinate, pointList[i - 1], pointList[i]) &&
                            IsOnLine(basel.Dest.Coordinate, pointList[i - 1], pointList[i]))
                        {
                            isget = true;
                            break;
                        }
                    }
                    if (!isget)
                    {
                        logger.Error("CheckIfError ConvexHull;basel:" + basel);
                        basel.LogOPrevEdge();
                        basel.Sym.LogONextEdge();
                        for (int i = 0; i < pointList.Length; i++)
                            logger.Info("ConvexHull:i:" + i + ";" + pointList[i]);
                        saveToProtobuf("D:\\abc", arr);
                        throw new Exception("CheckIfError basel not in ConvexHull");
                    }
                }
            }
            finally
            {
                //TimeComputeCL.TimeCompute.Pause("other");
            }
        }        

        private void saveToProtobuf(string filename, Coordinate[] contract)
        {
            try
            {
                using (var file = File.Create(filename))
                {
                    Serializer.Serialize(file, contract);
                }
            }
            catch(Exception e)
            {
                logger.Error("saveToProtobuf error:" + e.Message);
            }
        }

        private bool IsOnLine(Coordinate p, Coordinate orig,
            Coordinate dest)
        {
            decimal ccw = (dest.X - orig.X) * (p.Y - orig.Y) - (dest.Y - orig.Y) * (p.X - orig.X);
            if (ccw == 0 && p.X >= Math.Min(orig.X, dest.X) && p.X <= Math.Max(orig.X, dest.X)
                && p.Y >= Math.Min(orig.Y, dest.Y) && p.Y <= Math.Max(orig.Y, dest.Y))
                return true;
            else
                return false;
        }
        
        public static bool IsInCircumcircle(Coordinate A, Coordinate B, Coordinate C
            , Coordinate D)
        {
            decimal d2 = D.X * D.X + D.Y * D.Y;
            // aij 第i行，第j列
            decimal a11 = A.X - D.X;
            decimal a12 = A.Y - D.Y;
            decimal a13 = A.X * A.X + A.Y * A.Y - d2;
            decimal a21 = B.X - D.X;
            decimal a22 = B.Y - D.Y;
            decimal a23 = B.X * B.X + B.Y * B.Y - d2;
            decimal a31 = C.X - D.X;
            decimal a32 = C.Y - D.Y;
            decimal a33 = C.X * C.X + C.Y * C.Y - d2;
            decimal r = a11 * a22 * a33 + a12 * a23 * a31 + a13 * a21 * a32
                - a11 * a23 * a32 - a12 * a21 * a33 - a13 * a22 * a31;
            return r > 0M;
        }                

        public static Coordinate Circumcentre(Coordinate a, Coordinate b, Coordinate c)
        {
            decimal cx = c.X;
            decimal cy = c.Y;
            decimal ax = a.X - cx;
            decimal ay = a.Y - cy;
            decimal bx = b.X - cx;
            decimal by = b.Y - cy;

            decimal denom = 2 * Det(ax, ay, bx, by);
            if (denom == 0)
                throw new Exception("Circumcentre error;a:" + a + ";b:" + b + ";c:" + c);
            decimal numx = Det(ay, ax * ax + ay * ay, by, bx * bx + by * by);
            decimal numy = Det(ax, ax * ax + ay * ay, bx, bx * bx + by * by);

            decimal ccx = cx - numx / denom;
            decimal ccy = cy + numy / denom;

            return new Coordinate(ccx, ccy);
        }

        private double Det(double m00, double m01, double m10, double m11)
        {
            return m00 * m11 - m01 * m10;
        }

        private static decimal Det(decimal m00, decimal m01, decimal m10, decimal m11)
        {
            return m00 * m11 - m01 * m10;
        }        

        public static bool IsInCircleCC(Coordinate a, Coordinate b, Coordinate c,
                                        Coordinate p)
        {
            if (a.Equals(b) || b.Equals(c) || a.Equals(c) || p.Equals(a) || p.Equals(b) || p.Equals(c))
                return false;
            Coordinate cc = Circumcentre(a, b, c);
            decimal ccRadius = a.Distance(cc);
            decimal pRadiusDiff = p.Distance(cc) - ccRadius;
            return pRadiusDiff < 0M;
        }
    }
}
